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Introduction 

The propulsion element of the NASA Advanced Subsonic Technology (AST) initiative is 
directed towards increasing the overall efficiency of current aircraft engines. This effort 
requires an increase in the efficiency of various components, such as fans, compressors, 
turbines etc. Improvement in engine efficiency can be accomplished through the use of 
lighter materials, larger diameter fans and/or higher-pressure ratio compressors. However, 
each of these has the potential to result in aeroelastic problems such as flutter or forced 
response. To address the aeroelastic problems, the Structural Dynamics Branch of NASA 
Glenn has been involved in the development of numerical capabilities for analyzing the 
aeroelastic stability characteristics and forced response of wide chord fans, multi-stage 
compressors and turbines. 

In order to design an engine to safely perform a set of desired tasks, accurate information of 
the stresses on the blade during the entire cycle of blade motion is required. This requirement 
in turn demands that accurate knowledge of steady and unsteady blade loading is available. 
To obtain the steady and unsteady aerodynamic forces for the complex flows around the 
engine components, for the flow regimes encountered by the rotor, an advanced compressible 
Navier-Stokes solver is required. A finite volume based Navier-Stokes solver has been 
developed at Mississippi State University (MSU) for solving the flow field around multistage 
rotors. The focus of the current research effort, under NASA Cooperative Agreement NCC3- 
596 was on developing an aeroelastic analysis code (entitled TURBO-AE) based on the 
Navier-Stokes solver developed by MSU. The TURBO-AE code has been developed for 
flutter analysis of turbomachine components and delivered to NASA and its industry 
partners. The code has been verified, validated and is being applied by NASA Glenn and by 
aircraft engine manufacturers to analyze the aeroelastic stability characteristics of modem 
fans, compressors and turbines. 
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Summary of Accomplishments 


Version 4 of the TURBO- AE code was developed, checked and delivered to NASA and its 
industry partners. In this version of code, all possible inter-blade phase angles are analyzed 
using a single blade passage by incorporating phase-lagged boundary conditions. Both 
direct-store and Fourier-decomposed methods have been implemented into the code for the 
phase-lag analysis. In the direct-store method, all the relevant information are stored for 
lagging (time-shifting) the passage boundary conditions by the appropriate phase. Since the 
flow variables for the passage fluid boundaries need to be stored for the oscillation cycle, this 
method can become prohibitive in regards to memory requirements. Memory requirements 
on CRAY computers can be reduced by using the solid-state devices (ssds) to read and write 
instead of storing the variables within the core memory. This option, however, is not 
feasible within a work-station environment. To alleviate the problem on work-stations, a 
Fourier-decomposition analysis method was implemented. In this method, the time variation 
of the flow variables at the passage boundary was decomposed into their Fourier coefficients 
and only relevant coefficients were stored. This significantly reduces the storage 
requirement, however, the computational time is increased as Fourier decomposition as well 
as reconstruction of flow variables from the stored Fourier coefficients are required at each 
time step. 

Both of these methods were implemented into the TURBO-AE code. The code has been 
verified by applying it to a helical fan and released to NASA and its industry partners. Some 
of the results obtained from this version of the code are summarized in Refs. [1,2]. Further 
code validation results are summarized in Ref. [3]. 

Modification of the TURBO-AE, code in order to perform forced response calculations was 
started. The 3-D unsteady boundary conditions, developed by researchers at UTRC [2], 
currently does not allow disturbances to enter the computational domain. In other words, it 
only allows disturbances to propagate out of the domain. Work was also started with the 
ultimate aim of understanding the unsteady boundary conditions and modifying them for 
forced response calculations. This modification requires coding changes to allow for 
prescribed disturbances from outside the domain to propagate into the computational 
domain. 

Drs. R. Srivastava and M. A. Bakhle, both of the University of Toledo, held two workshops, 
for the industry partners and the Air Force, at NASA Glenn. The workshops provided the 
potential users of TURBO-AE all the relevant information required to prepare the input data, 
execute the code, interpret the results and benchmark the code on their computer systems. 
Technical support was also provided to researchers at NASA and to the industry partners 
who are currently using the code. 
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ABSTRACT 

In the present work a comparative study of phase- 
lagged boundary condition methods is carried out. 
The relative merits and advantages of time-shifted 
and the Fourier decomposition methods are compared. 
Both methods are implemented in a time marching 
Euler /Navier-Stokes solver and are applied to a flat plate 
helical fan with harmonically oscillating blades to perform 
the study. Results were obtained for subsonic as well as 
supersonic inflows. Results for subsonic inflow showed 
good comparisons with published results and between the 
two methods along with comparable computational costs. 
For the supersonic inflow, despite the presence of shocks 
at the periodic boundary results from both the methods 
compared well, however, Fourier decomposition method 
was computationally more expensive. For linear fiowfield 
Fourier decomposition method is best suited, especially 
for work-station environment. The time-shifted method 
is better suited for CRAY category of computers where 
fast input-output devices are available. 

INTRODUCTION 

Numerical aeroelastic analysis methods have been 
developed using both two-dimensional and three- 
dimensional aerodynamics. Srivastava et al. (1998) have 
provided a good reference for these methods. The meth- 
ods based on two-dimensional aerodynamics are fast but 
ignore the real physical effects of three-dimensional flows. 
The three-dimensional analyses capture all the required 
physics but are much more computationally expensive. 
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More so for analyses of inter-blade phase angles (IBPAs) 
requiring large number of blade passages to be included 
in the analysis. To reduce computational time, especially 
for smaller IBPA vibrations, several phase-lagged meth- 
ods have been reported in the past. Erdos et al. (1977) 
were the first to develop an analysis based on direct- 
store method. The direct-store method stores all the 
relevant fluid properties over the oscillation cycle which 
is applied with appropriate lag for the IBPA being an- 
alyzed. Though, no loss in fidelity occurs, the method 
requires significant additional memory and becomes pro- 
hibitive for large three-dimensional problems. Another 
method, known as time-inclined computational plane ap- 
proach, was proposed by Giles (1983) primarily to over- 
come the problem encountered in rotor-stator applica- 
tions where no final periodic state exists. This method 
requires transforming the original governing equation to 
account for the tilting of the time plane. He (1989) pro- 
posed a shape-correction method for applying the phase- 
lagged boundary conditions that did not have the storage 
penalty associated with the direct-store method. In the 
shape-correction method, the variation of fluid properties 
over an oscillation cycle is decomposed into its Fourier 
coefficients and only the coefficients are stored. These co- 
efficients are used later to regenerate the fluid properties 
as required. Later, He and Denton (1994) extended the 
method to three-dimensions. Peitsch, Gailus and Weber 
(1996) proposed a variation of the direct-store method 
to reduce the storage requirements, using a “foothold 
technique” that stores the fluid properties only at cer- 
tain foothold points. The properties at other points over 
the oscillation cycle is obtained by interpolation from the 
nearest foothold point. 

Except for the Giles’ “time-tilting” method, which re- 
quires transforming the governing equations, the above 
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methods are based on either the direct-store method or 
the shape-correction method. The direct-store method 
requires large in core memory for storing the fiowfield 
properties over the oscillation cycle. On a CRAY type 
machine where solid state devices (SSDs) can be used for 
fast input-output, the memory requirements of the direct- 
store method can be minimized. The Fourier decompo- 
sition method does not require large data storage, how- 
ever, additional computational time is required to Fourier 
decompose and then regenerate the fluid properties us- 
ing stored coefficients. Unlike the direct-store method, 
it relies on superposition making the problem essentially 
linear. Linearity may be of concern for problems where 
strong vibrating shocks are present at the periodic pas- 
sage boundary. Further, because of the lag associated in 
enforcing the ‘‘phas e-lag” , the rate of convergence also 
becomes an issue for the two methods. Clearly, there are 
advantages and disadvantages associated with both the 
methods. It is not clear if one method is superior to the 
other. The present study attempts to highlight the advan- 
tages of one method over the other, given the problem of 
interest and resources at hand. Towards this goal, the pri- 
mary objective of the present study is to implement both 
phase-lagged methods into one solver. Both the methods 
are then applied to identical problems in order to investi- 
gate their relative advantages. 

In the present work, the two methods are implemented 
within the TURBO-AE code. The TURBO-AE code is 
currently under development at NASA Lewis Research 
Center. The details of the code along with several results 
have been reported by the authors of this paper in Bakhle 
et al (1996, 1997). This analysis can analyze flutter for all 
the possible IBPAs. Srivastava et al. (1998) implemented 
the time-shifted method based on direct-store method and 
validated and verified the analysis with previously pub- 
lished results. In the present work the Fourier decompo- 
sition method with multiple updates per os cill ation cycle 
is implemented within the TURBO-AE and the progr am 
is applied to a flat plate helical fan. The obtained results 
are compared for accuracy and computational efficiency 
for each of the two methods. To reduce computational 
cost, results are obtained using invisrid calculations. 

THE TURBO-AE CODE 

The aeroelastic solver TURBO-AE is described in brief 
in this section. Interested readers may refer to Bakhle 
et al (1996, 1997) for greater details. The details of the 
time-shifted boundary conditions are described in Srivas- 
tava et al. (1998). The Fourier decomposition method is 
based on the method proposed by He (1989). TURBO- 
AE is based on an Euler /Navier-S tokes unsteady aero- 
dynamic solver TURBO (Janus 1989), for internal flow 


calculations of axial flow turbomachinery' components. 
TURBO-AE can model multiple blade rows undergoing 
harmonic oscillations with arbitrary IBPAs. The aerody- 
namic loads are obtained by solving the unsteady Eu- 
ler or Navier-S tokes equations. For the viscous calcu- 
lations, Reynolds-averaged Navier-Stokes equations are 
solved. The Baldwin-Lomax equations are used to model 
the turbulence. The aerodynamic equations are solved us- 
ing a finite volume scheme. Flux vector splitting is used 
to evaluate the flux Jacobians on the left hand side. The 
right hand side fluxes are discretized using the higher or- 
der Toted Variation Diminishing (TVD) scheme based on 
Roe’s flux difference splitting. Newton sub-iterations are 
used at each time step to maintain the higher accuracy. 
Symmetric Gauss-Sidel iterations are applied to the dis- 
cretized equations. A newly developed three-dimensional 
non reflecting boundary condition (Montgomery and Ver- 
don 1997) is applied at the upstream and downstream 
boundary. The blade motions are at a prescribed fre- 
quency, and are simulated using a dynamic grid deforma- 
tion technique. The grid is updated at each time step by 
recalculating the grid using linear interpolation, assuming 
the far field boundaries to be fixed. The grids on the cas- 
ing are, however, allowed to slide along the casing. The 
aeroelastic characteristics of the rotor are obtained by cal- 
culating the energy exchange between the vibrating blade 
and its surrounding fluid. A positive work on the blade 
indicates instability (Bakhle et al. 1996). 

NUMERICAL RESULTS 

Sample results obtained from TURBO-AE for phase- 
lagged boundary conditions axe presented in this section. 
Results are obtained for a flat plate helical fan config- 
uration used by Montgomery and Verdon (1997). The 
fan configuration consists of 24 flat plate blades with zero 
thickness enclosed within a rigid cylindrical duct with no 
tip-gap. At mid-span, the stagger angle is 45° and the gap 
to chord ratio is one. Results axe presented for two dif- 
ferent inflow conditions at mid- span: a subsonic relative 
inflow Mach number of 0.7 at zero incidence with axial 
Mach number of 0.495 and a supersonic relative inflow 
Mach number of 1.3 at zero incidence with axial Mach 
number of 0.9192. 

The grid used for the analysis is an H-0 type grid with 
141 points in the streamwise direction, 11 points in the 
spanwise direction and 41 points in the blade to blade 
direction. The aeroelastic analysis is carried out by first 
obtaining a steady aerodynamic solution for the given con- 
ditions. From this steady solution, the unsteady solution 
is started by forcing the blades to undergo a harmonic mo- 
tion at the given frequency, mode shape and IB PA. The 
unsteady aerodynamic behavior, as well as work-per-cycle 
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unsteady pressure difference 


is calculated for the oscillating blades. 

A comparison of results from TURBO- AE, a linearized 
Euler analysis (Montgomery and Verdon 1997) and a two 
dimensional linear theory (Smith 1972) is shown in Figs. 
1 <k 2 at the subsonic relative inflow condition. These 
results are reproduced here from Srivastava at al (1998) 
for sake of completeness. A good comparison with lin- 
ear theory and linearized analysis indicates the unsteady 
aerodynamic behavior predicted by TURBO- AE to be ac- 
curate. 



Figure 1: Unsteady pressure difference variation with 
chord at mid-span for 0 deg IBPA pitching oscillations 

The phase-lagged analysis was carried out next for the 
subsonic relative inflow condition. The blades were forced 
into a pitching oscillation about their mid-chord at a re- 
duced frequency of one and -90 deg IBPA. The analy- 
sis was earned out using both the time-shifted and the 
Fourier decomposition methods. Two different analyses 
were performed for the Fourier decomposition method. 
In one of the analyses one Fourier coefficient was retained 
and the coefficient was updated only once per oscillation 
cycle. In the other analysis the coefficient was updated 
four times during each oscillation cycle. The multiple up- 
dates of the coefficients is expected to provide a faster 
convergence. The time-shifted analysis has been previ- 
ously validated by Srivastava et al (1998). The results 
obtained for the Fourier decomposition method are, there- 
fore, compared with the results obtained from the time- 
shifted method to ascertain the accuracy of the Fourier 
decomposition method. These results are shown in Figs. 
3 & 4. In Fig. 3 the variation of total work-per-cycle with 
oscillation cycle is shown. It can be seen that the Fourier- 



Figure 2: Unsteady pressure difference variation with 
chord at mid-span for 180 deg IBPA plunging oscillations 


decomposition method with multiple updates per cycle is 
the fastest to converge with single update per cycle be- 
ing the slowest. However, the three analyses eventually 
converge within 0.1 % of each other, indicating that once 
convergence is achieved the results are same. This is fur- 
ther verified by comparing the unsteady pressure differ- 
ence variation along chord at mid-span. The comparison 
for the three analyses is shown in Fig. 4. A very good 
comparison is obtained. This indicates that the results 
from the Fourier decomposition method are as accurate 
as the time-shifted analysis. It also indicates that the 
four updates per cycle converges much faster than single 
update, hence in all subsequent work four updates per os- 
cillation cycle will be used for the Fourier-decomposition 
method. 

For the subsonic inflow condition, the analysis was also 
carried out using four blade-passages in order to simulate 
the -90° IBPA condition. Using four blade-passage anal- 
ysis provides the most accurate results as no approxima- 
tions are involved. The phase lagged conditions introduce 
errors into the analysis especially since the conditions at 
the start of the analysis are not known. Comparisons with 
four-blade-passage analysis also provide a means to mea- 
sure the benefits of phase-lagged boundary conditions in 
terms of savings of computer resources. The results ob- 
tained for the four-passage analysis are compared with the 
results obtained from the time-shifted and the Fourier- 
decomposition methods. The convergence of work-per- 
cycle for the three analysis methods is shown in Fig. 5. 
For the four passage analysis, as expected, the total work 
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Figure 3: Comparison of work-per-cycle convergence for 
phase-lagged analysis methods for M„=0.7 and -90 deg 



Figure 5: Comparison of work-per-cycle convergence for 
phase-lagged analyses with multiple passage analysis for 
Afoo =0.7 and -90 deg IB PA 



Figure 4: Comparison of unsteady pressure difference 
variation with chord at mid-span for phase-lagged analy- 
ses for Moo =0.7 and -90 deg IBPA 


for all the four passages coalesce. Also, the four passage 
analysis takes approximately four to five oscillation cy- 
cles to converge. The time-shifted analysis, shown with 
dashed lines, requires approximately eight cycles, whereas 
the Fourier-decomposition method converged in five to six 
oscillation cycles. The total work from the three analysis 
methods converge to within 0.5% of each other. This in- 
dicates that the three methods are equally accurate. This 
is. further confirmed by comparing the unsteady pressure 
difference variation, Fig. 6. Once again a very good com- 
parison is obtained indicating accuracy is not a concern 
for phase-lagged boundary conditions for the conditions 
where flowfield can be assumed linear. 

FVom these results it can be seen that the four-passage 
analysis requires the least number of oscillation cycles to 
converge, whereas the time-shifted analysis requires the 
largest number of cycles. However, the computational 
cost for the four-passage analysis was largest. This is 
because the analysis had to be carried out using four- 
passages as opposed to a single passage for phase-lagged 
analyses. This reduced the problem size of the phase- 
lagged analysis to one fourth that of the four-passage 
analysis. The smaller number of cycles required for con- 
vergence do not sufficiently offset the increase in com- 
putational cost for the multiple passage analysis. Fur- 
ther, despite the difference in rate of convergence for 
the time-shifted analysis and the Fourier-decomposition 
method, the computational costs required for convergence 
are fairly comparable, see Table 1. This is because of the 
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unsteady pressure difference 



Figure 6: Comparison of unsteady pressure difference 
variation with chord at mid-span of phase-lagged anal- 
yses with multiple passage analysis for M 00 = 0.7 and -90 
deg IB PA 


increased computational cost required by the Fourier de- 
composition method to Founer decompose and regenerate 
the fluid properties is offset by increased rate of conver- 
gence. It should also be noted here that the memory re- 
quirements of the time-shifted analysis can be reduced by 
using the SSDs available on the CRAY computers. The 
SSDs help reduce the memory requirements sig n ifi can tly 
for a marginal increase in the computation cost of reading 
and writing to the disk. 

All of the above analyses were carried out on a CRAY 
C-90 computer. The computer resources required are tab- 
ulated in Table 1. The CPU time required per time step 
shows the additional cost per time-step for the Fourier de- 
composition method over the time-shifted method. This 
increase in time, however, is more than compensated for 
by increased rate of convergence. Also, for the current 
problem, using the Fourier decomposition method reduces 
the memory by over 50% as compared to the time-shifted 
analysis. Also, it can be seen that despite the faster rate 
of convergence for the four-passage analysis the CPU re- 
quirements are almost three times as much as that of the 
phase-lagged methods. This difference will increase sig- 
nificantly for analyzing the smaller IBP As requiring many 
more blade passages for the multiple passage analysis. Es- 
pecially since the rate of convergence for the time-shifted 
analysis is independent of the IBPA being analyzed (Sri- 
vastava tt oi, 1998). Also, from Table 1 it appears that the 
Fourier-decomposition method with multiple updates of 


the coefficient and the time-shifted method may be com- 
parable in CPU requirements. Therefore, the choice of the 
particular method will depend on the available resources 
and the problem being analyzed. For problems where flow 
behavior can be assumed linear, superposition is not of 
concern, the Fourier-decomposition method is best suited 
irrespective of the resources available. For problems where 
nonlinearity may be of concern time-shifted method may 
have to be used, especially if SSDs are available. 

The two methods were next applied to a supersonic in- 
flow condition to help evaluate the effectiveness and prob- 
lems that Fourier-decomposition method might have for 
flows with nonlinearities. The analysis was carried out for 
+90 deg and -90 deg IBPA. The flow condition is subres- 
onant (waves propagating for supersonic relative inflow, 
Verdon (1989)) for the +90 deg IBPA and is superres- 
onant (waves decay away from the blade row, Verdon 
(1989)) for the -90 deg IBPA with resonance at -102.1 
deg. The +90 deg was analyzed using the Fourier de- 
composition method with one, five, and 11 Fourier coef- 
ficients, as well as with time-shifted method. A Fourier- 
decomposition analysis was carried out using 15 coeffi- 
cients also, but the results were found to be identical to 
the 11 coefficient analysis. 

The comparison of the total work convergence history 
for the four analyses is shown in Fig. 7. For sake of com- 
putational cost the analysis was carried out for only 18 
oscillation cycles. Even though the flowfield is not com- 
pletely converged after 18 cycles, the disturbances appear 
to be dying out. The 11 coefficient analysis compares very 
well with the time-shifted analysis. 

The one and five coefficient analyses show small differ- 
ences. Interestingly the total work calculated using only 
one coefficient compares very favorably with time-shifted 
analysis and also converges faster. However, significant 
differences are found for the unsteady pressures between 
the results from one coefficient analysis and the other 
three analyses. The first and second harmonics of the 
blade surface unsteady pressure differences at mid span 
are shown in Fig. 8. The first harmonic pressure shows 
good comparison between the four analyses over most of 
the blade chord. Over the last 10-12 % of the chord, 
one coefficient analysis shows differences in both real and 
imaginary pressures. The second harmonic of the pressure 
difference, on the other hand, shows significant differences 
between the one coefficient Fourier analysis and the other 
three analyses. The time-shifted, five and 11 coefficient 
analyses show good comparison with each other over most 
of the chord, with some minor differences near the trail- 
ing edge for the five coefficient analysis. The 11 coefficient 
analysis compares very well with the time-shifted method, 
hence in all future analyses 11 coefficients were used. 
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The analysis was next carried out for the -90 deg IBP; 

tTmlTf T °\ ^ 15 Sh0WD b Fi S- 9 fo: r hot 

tune-shifted analysis and the Fourier decomposition ana 

ysis using 11 coefficients. The analysis was again cai 

ned out for 18 vibration cycles. The time-shifted metho 

convergence, whereas the Fourier decompositio 

rion ^ d p 0eS n °i‘ ThUS that for this condi 

tion the Founer decomposition method is computational! 

more expensive than the time-shifted method. The tota 
work from both the methods are comparable even thoug] 
decom P°sition analysis has not totally con 
verge . The first and second harmonics of the unstead 1 

pressure difference at mid span are shown in Fig 10 Th< 

companson between the two methods is good except ova 
the last 10-15% of the blade chord near the trailing edge 
During parts of the oscillation cycle a shock appears in 
this region of the blade along with passage shocks at the 
uid periodic boundaries upstream of the blades. The 
owfield also showed the passage to be choked for parts 

ThCSe fl ° W featur es indicate that the 
^wfidd behimor may be nonlinear. For the +90 deg 
1BPA, shocks were present but neither choked flow nor 
passage shocks at the fluid periodic boundaries were ob- 
served. Despite the presence of shock and choked flow 
condxtions the method of superposition provides compa- 
rabk results for this configuration. However, it should be 
noted here that the geometry without loading and flow 
urning is simplistic. Further investigations with realistic 
geometry needs to be carried out to understand in more 
detail the influence of nonlinearities on the Fourier de- 
composition method. 


(a) First Harmonic 



(b) Second Harmonic 


Figure 8: Comparison of unsteady pressure difference 
variation with chord for M^l.3 and +90 deg IBPA 
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Figure 9: Comparison of work-per-cycle history for -90 
deg D3PA at M^l.3 and -90 deg EBPA 

(a) First Harmonic 


CONCLUSIONS 

Two phase-lagged boundary condition methods have 
been successfully implemented into the TURBO-AE, an 
Euler /Navi er- Stokes based aeroelastic solver. Both these 
methods along with the multiple passage analysis method 
have been applied to an identical geometry to investigate 
the accuracy and efficiency of the various methods. Com- 
paring the results obtained from these methods it was 
found that all the methods provide equally accurate re- 
sults, for the subsonic relative inflow condition. For the 
supersonic relative inflow the presence of a shock was in 
itself not sufficient to invalidate the Fourier decomposi- 
tion method. The errors introduced by superposition were 
negligible if large number of Fourier coefficients were in- 
cluded in the analysis even for flows with shocks at the 
periodic boundaries and choked flow conditions. 

The study also showed that the phase lagged methods 
significantly reduce the computational cost as compared 
to the multiple passage analysis. These savings will be 
much more significant for smaller EBPAs. The CPU cost 
of the time-shifted method was found to be comparable 
to that of the Fourier decomposition method. However, 
because of the large reduction in memory requirements, 
it is recommended that for flowfields where superposition 
is not of concern, Fourier-decomposition method be used. 
For flows with nonlinearities more study is required to 
further quantify the nature of the flowfield for which the 
Fourier decomposition method may break down. 



(b) Second Harmonic 


Figure 10: Comparison of unsteady pressure difference 
variation with chord for Moo =1.3 and -90 deg IBPA 
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Version 

CPU time/ 
time-step 

CPU time for 
convergence 

CPU Mem. 

4 Passages, in core 
storage 

27.36 secs 

7 hrs 45 mins 

67 Mws 

Time-Shifted, in core 
storage 

6.94 secs 

L 

2 hrs 30 mins 

34 Mws 

Time-Shifted, ssds 
used for I/O of BCs 

6.98 secs 

2 hrs 35 mins 

15 Mws 

Fourier-Decomposition 
1 Coefficent, 1 update 

7.29 secs 

4 hrs 

15 Mws 

Fourier-Decomposition 
1 Coefficent, 4 updates 

7.67 secs 

2 hrs 15 mins 

16 Mws 


Table 1: Comparison of computer resources required by various methods for Af co =0.7 and -90 deg IBPA 
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ABSTRACT 

In the present work the unsteady aerodynamic char- 
acteristics of harmonically oscillating fan blades are in- 
vestigated by applying a time-shifted boundary condition 
at the periodic boundaries. The direct-store method is 
used to implement the time-shifted boundary condition 
in a time-marching Euler/Navier-Stokes solver. Invisdd 
flow calculations for a flat plate helical fan, in a single- 
blade passage domain, are used to verify the analysis. 
The results obtained show good correlation with other 
published results as well as with the same solver using 
multiple blade passages stacked together. Significant sav- 
ings in computer time is realized, especially for smaller 
phase angles. 

INTRODUCTION 

The objectives of the Advanced Subsonic Technology 
(AST) project, funded by NASA, is to improve the effi- 
ciency of the turbomachines and reduce the NOX emis- 
sions. To satisfy these objectives, numerical aeroelastic 
analysis methods for turbomachinery applications are cur- 
rently being developed at NASA Lewis Research Cen- 
ter. Numerical aeroelastic analysis methods have pri- 
marily been developed for two-dimensional cascades, rep- 
resentative of turbomachinery configuration (for exam- 
ple, Platzer 1978, Sisto 1977, Fleeter 1979, and Bendik- 
sen 1990, among others). These methods are inade- 
quate for an accurate analysis as they ignore the strong 
three-dimensional flow characteristics present in a turbo- 
machine. The large computational cost associated with 

•Senior Research Associate, also Resident Research Associate, 
NASA Lewis Research Center 
Distinguished University Professor 
Machine Dynamics Branch 


three-dimensional analysis, has restricted the primary fo- 
cus for the three-dimensional aeroelastic analyses, over 
the years, to predicting unsteady aerodynamic forces, 
assuming the unsteady behavior to be linear. Several' 
methods for calculating the unsteady aerodynamic loads 
over vibrating blades have been reported in the literature 
by solving linear potential equations (Namba 1987, Chi 
1993), linearized Euler equations (Hall and Lorence 1993, 
Montgomery and Verdon 1997) and non-linear Euler and 
Navier-Stokes equations (He and Denton 1994, Peitsch, 
Gallus and Weber 1996, and Bakhle et a l. 1996 k 1997). 
A good review of analytical methods for turbo machin ery 
blade vibrations is given by Chi (1993). 

Some three-dimensional aeroelastic analyses of turbo- 
machinery components have also been reported. Carta 
(1967), used a quasi three-dimensional approach by stack- 
ing two-dimensional strips of isolated airfoil. W illiam* & 
Cho (1991) reported an analysis based on linear panel 
method and solved the aeroelastic equation using a pulse 
and influence coefficient approach. Recently, methods 
based on Euler analysis have been reported by Geroly- 
mos (1992), Srivastava and Reddy (1995) and Srivastava, 
Reddy and Stefko (1996). Srivastava and Reddy (1997) 
have also applied various aeroelastic analysis t echni ques 
for flutter analysis of a ducted rotor configuration and in- 
vestigated the advantages and disadvantages of <»a r h of 
the methods. 

With the exception of He and Denton (1994) and 
Bakhle et al (1997), all the other methods reported in 
the literature, ignore the effects of viscosity. He and Den- 
ton (1994) have used thin layer approximation to include 
the effects of viscosity in their calculations. To accu- 
rately model flows with separation due to stall or shock 
boundary layer interaction an aeroelastic analysis pro- 
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dased 011 full Navier-Stokes equations, 
currently being developed at NASA Lewis Research 
enter. The numerical analysis within the TURBO- AE 
program, couples an aerodynamic analysis based on three- 
dimensiona 1 unsteady Euler/Navier-Stokes equations with 
anormal mode structural analysis. The aeroelastic char- 
acteristics are obtained by calculating the energy ex- 
change between the rotor and the surrounding fluid. 

er JrLuft TURB °- AE P r °S™m along with sev- 

al results have been reported by the authors of this pa- 
per m Bakhle et al (1996, 1997). This analysis can L 

tmWTh f ° r f the P° ssible **«■*** P^e angle 
(IBPA). The number of IBPAs possible for any given n> 

1956) Sr ' t ^ thC ““t" ° f blad6S “ the rot ° r ( L ane 
sS the tben ° D - zero L IBp A calculations, the analysis 

tSto t nUmber ° f blade t0 simu late 

he motion with appropriate IBPA. This is a cumbersome 

StheSnl C f 0nSUmm S P rocess - To improve the efficiency 
L th ^°lT\ 3 tune ' s hifted boundary condition is bZ 
mg added to the TURBO-AE solver. In this method all 

£TbS Ie BPA "to**™ be Performed XI 

srngle blade passage by applying a time-shifted boundary 
condition across the periodic boundaries of the passag^ 

Several researchers have reported methods in two- 
dunensions for using time-shifted periodic boundary con- 

1197? r l UC r he C ° mputational domain. Erdos^ et al. 
(1977) were the first to develop a method based on direct 

oftZZtZf- ^/‘-ative approach to the protS 

by XX^The" C0 ? m ° n developed 

aDnrn . , ' , J/' The tlme - inclined computational plane 
P ° Gdes was Primarily to overcome the problem 
encountered m rotor-stator applications, where no final 

s “ e ^ He < i989 > ■ <» °'<*« zz.* 

storage requirement of the direct-store method, proposed 

X r O fT tl0n meth ° di Wherein ’ onl y Fourier co 
at th! ? ° f J e unstead - V variation of the fluid properties 
fl9«q's Pm0dlC b ° Undar >’ were stored. The method of He 
ton f 19941 e ? e “ d t d fc ° three * dimensions by He and Den- 
r eD oL d 4 ^ Pei £ s< J> GaU us and Weber (1996) have also 
for Sthrl 6 time-shifted boundary condition 

biSTn H ,^ enS1 ° nal Euler S0lver - Their method is 
direct-store method, but in order to reduce the 

t^T®? 8 a “ footbold technique” is used. L 
th T* f ^ the fluid properties over the entire cycle 

points°The ieS are r ! t0red at a smaU number of foothold 
pomts. The properties at other time steps during the cv- 

poims ° tamCd ^ mter P olat ion from the nearest foothold 


tio^h^r 63 ^ W ° rk ’ tbC time-shifted boundary condi- 
on has been incorporated within the TURBO-AE nro- 
gam usmg the direct-store method. To verify aTd Z 
he solver, it is apphed to a flat plate helical fan 


operating in subsonic conditions. This provides for a res 
sonably good test case as there are some numerical result 
that have been published for this geometry. The flat plat 
helical fan, and the operating conditions are such that th 
mid span section of the fan has very near two-dimensiona 
flow, at zero incidence, over it. Results obtained for th, 
fan are compared with numerical results from linear the 
ory (Smith 1972) and linearized Euler analysis (Mont 
gomery and Verdon 1997), as well as with resuifob^ 

TjZTr^T ° PtiOQ of ' TURBO-AE (Bakhlt 
et al 1997). To reduce the computational cost, the results 
are obtained usmg mvisdd calculations. 


THE TURBO-AE CODE 

. The aero das tic solver TURBO-AE is described in bri, 

“ ?no™ 10 ?™!? tereSted readers may refer to Bakhl 
(1996, 199 f) for greater details. TURBO-AE i 

baSed “ Euler/Navier-Stokes unsteady aerodynami 
solver TURBO (Janus 1989), for internal flow ikufe 
ions of axial flow turbomachinery components. TURBO 

osdll^? model f lulti P le blade rows undergoing harmonii 
escalations with arbitrary IBPAs. The aerodynami, 

loads are obtained by solving the unsteady Euler o, 

R^lf Rations. For the viscous calculations, 
Reynolds-averaged Navier-Stokes equations are solved, 
ihe Baldwm-Lomax equations are used to model the tur- 
ence. The aerodynamic equations are solved using a 
finite volume scheme. Flux vector splitting is used to 
evaluate the flux Jacobians on the left hand side. The 
nghthandsKie fluxes are discretized using the higher or- 
der Total Vanation Diminishing (TVD) scheme based on 
Roe s flux ffifference sphtting. Newton sub-iterations are 
used at each time step to maintain the higher accuracy. 
Symmetric Gauss-Sidel iterations are applied to the dis- 
cretized equations. A newly developed three-dimensional 
non reflecting boundary condition (Montgomery and Ver 
don 1997) is applied at the upstream id 2^ 
boundary. The blade motions are at a prescribed fre- 
quency, and are simulated using dynamic grid deforma- 
tion technique. The grid is updated at each time step by 
recalculating the grid using linear interpolation, ass uming 
the far field boundaries to be fixed. The grids on the cas- 
ing are, however, allowed to slide along the casing. The 
aeroelastic characteristics of the rotor are obtained by cal- 
culating the energy exchange between the vibrating blade 
and its surrounding fluid. A positive work on the blade 
indicates instability. 


TIME-SHIFTED BOUNDARY CONDITIONS 
In the work reported by Bakhle et al (1997), the non 
zero IBPA analysis was carried out by stacking the re- 
quired. number of blade passages. For smaller IBPA anal- 
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ysis, depending upon the number of blades in the rotor 
the required number of passages could be very large and 
computationally prohibitive. To alleviate this problem 
a time-shifted boundary condition based on direct-store 
method, has been implemented into the analysis. Since 
the time-shifted boundary condition analysis is carried 
out using a single passage for any IBPA, the computa- 
tional cost is significantly reduced, especially for smaller 
JLdFA analysis. 

. ^ he f ourier decomposition of the fluid properties, as 
m the shape-correction method of He and Denton ( 1994 ) 
assumes linearity of the fluid properties. This could lead 
to inaccuracies for strongly non-linear flows at the peri- 
odic bomidary such as presence of a shock. Also, since 
the coefficients are calculated at the end of each cycle 

the application of the boundary conditions are lagged by 
one full cycle. 66 y 

In the present analysis, the periodic boundary condi- 
tions are lagged only by the phase on one surface and by 
the difference of 360° and phase on the other. The reduced 
lag in boundary condition application helps achieve con- 

At the be S inn “g of &e analysis, where 
the lagged fluid properties are not available, the instanta- 
neous values are used. This introduces some initial errors 
that are driven out of the computational domain within 
a few cycles of oscillation. These initial errors result in 
increased number of oscillations for convergence as com- 
pared to the original analysis with multiple blade pas- 
sages. The large storage requirements of the direct-store 
method can be mi nimiz ed by writing the fluid properties 
on to the sohd state devices (ssds), rather than storing 
them in core memory. The ssds on CRAY C-90 allow for 
large input-outputs without any appreciable degradation 
in code performance. 


NUMERICAL RESULTS 
Sample results obtained from TURBO- AE time-shifted 
anatysis for the flat plate helical fan configuration, used 
by Montgomery and Verdon (1997), are presented in this 
section. The helical fan consists of flat plate blades with 
zero thickness enclosed within a rigid cylindrical duct with 
no gap between the blade tip and the duct. The blades are 
tested to maintain a zero incidence at all sections for the 
inflow. The inflow relative Mach number at midspan is 
0.7 with the free stream axial Mach number being 0.495. 

he stagger angle at mid-span of the cascade is 45°, and 
the gap-to-chord ratio is unity for a rotor with 24 blades 
and a diameter of 8.448. The hub to tip ratio is 0.8. 

i ^ he - grid . USed for the ^^ysis is an H-0 type grid with 
141 pomts m the streamwise direction, 11 grid points in 
the span wise direction and 41 points in the blade to blade 
direction. The aeroelastic analysis is carried out by first 


obtaining a steady aerodynamic solution for the given con 
ditions. From this steady solution, the unsteady solution 
is started by forcing the blades to undergo the harmonic 
motion at the given frequency, mode shape and IBPA. 
The unsteady aerodynamic behavior, as well as work-per- 
cycle is calculated for the oscillating blades. 

In the first few calculations, the code is verified by ob- 
taining the unsteady behavior of the blades undergoing ei- 
ther pure pitching motion or pure plunging motion The 
plunging motion is normal to the blade chord at mid 
span and is of constant amplitude for the entire span 
The pitching motion is about the mid chord. The non’ 
dmaensional frequency of oscillation, based on blade chord 
for these analyses has been taken as unity. An oscilla- 
tion amplitude of 0.2° for pitching and 0.1% of chord for 
plunging was used. It was found that 200 steps per os- 
ailation cycle were required to eliminate dependency on 
tune steps. Also, for the multiple passage analysis, it was 
lound that a minimum of four to five oscillation cycles 
were required to obtain a converged solution. The work 
per-cycle convergence is shown in Fig. 1 for three different 
time steps. Reducing the time step from 100 steps per os- 
cillation cycle to 200 steps, a large difference in work is 
seen. However, reducing the time step further does not 
impact the solution significantly. Therefore, in all subse- 
quent calculations, 200 time steps per oscillation cycle are 
used. 

Figures 2 - 5 show the variation of the unsteady pressure 
along the chord at mid-span. The variation of the first 
harmonic of real and imaginary parts of the difference of 
the unsteady pressure between the pressure and suction 
surfaces is plotted against the chord. These results are 
obtained without using the time-shifted boundary condi- 
tions and are presented here for validation purposes of the 
imsteady calculations of the code. For the non-zero IBPA. 
the required number of passages were stacked for these 

“►wf \ ^ be S6€n ’ good COI *Parison is obtained 

withthelmeanzed Euler analysis (Montgomery and Ver- 

on 1997) and the analytical results (Smith 1972). These 
figures indicate that the unsteady behavior of the solver 
compares well with other published results. 

To obtain the time-shifted results, the analysis was 
earned out for various IBPAs for both the pitching and 
plunging cases. Some of these results are presented here 
for verifying and validating the method. The total work at 
the end of each cycle was monitored to obtain the conver- 
gence. The Ration of total work-per-cycle with vibra- 
lon cycle for -90° IBPA is shown in Fig. 6. Also, shown in 
this figure is the variation obtained from the four passage 
analysis. It can be seen that the time-shifted analysis re- 
quires several more vibration cycles to reach convergence. 
Typically, the multiple passage analysis showed conver- 
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gencern 4-o cycles, whereas, the time-shifted analysis re- 
quired ' -10 cycles to reach convergence. This is to be 
expected. For the multiple passage analysis the bound- 
b ? m the start applied appropriately 
M h ^ flUld , mterfaces ' 0n the o^er hand, for the time- 
shifted analysis, the lagged fluid properties at the begin- 
ning of the solution are not available, hence, inaccurate 
boundary condition is applied at the start of the analysis. 
The errors introduced because of this inherent drawback 
ot the method, require longer to move out of the calcu- 
lation domain, thus requiring more cycles of oscillations 
for convergence. From this figure one can also see that 
the work obtained from all the passages of the multiple 
passage analysis are same once the solution has reached 
convergence. It should be noted here that, even though, 
e m tiple passage analysis requires approximately half 
e oscillation cycles for convergence, it has to use four 
lade passages. This in turn results in the multiple pas- 

sageanalysis requiring approximately twice the CPU time 
tor this case. 

The real and imaginary parts of the unsteady pressure 
difference for the two methods are compared with each 
other in Figs. 7 & 8 for 180° IBPA for pitching and plung- 
ing motion, and in Figs. 9 & 10 for 90° and -90° IBPAs for 
pitching oscillation. Also shown in Figs. 9 & 10 are results 
from linearized Euler (Montgomery and Verdon 1977) and 
linear theory (Smith 1972) for comparison purposes. In 

mLtf 2®“ f° d com P arison i 5 obtained, in fact for 

f“u,° m' reSUltS *" 

problem was encountered in conver- 
gence. This case is a superresonant condition and hence 
had some reflections from the upstream boundary. Even 
though non-reflecting boundary conditions are used, the 
radial modes are not accounted for in the analysis. This 

“ SOme ^Sections from the boundary. Because 
of this the computations required longer for convergence. 
Tins was observed for both the multiple passage analysis 
well as the tune-shifted boundary condition case. 

The analysis was carried out for 45° and -45° IBPA 
as well. These results also compared well with the lin- 

tTJn \° f Smith ( - 1972) - ° De “ teres ting fact seems 
Of F w rge b T.i°“ panng the conv ergence characteristics 
the mpT e 'f hlfted b0Undary ^o'dations - irrespective of 
ffiPA of motion, the convergence is achieved within 
seven to 10 oscillation cycles. This is shown in Fig U 
This unphes that potential savings are quite large for 
smaller IBPAs since they require larger number of pas- 
sages for the multiple passage analysis. The break even 

mpV tu enns ° f com P utati °nal costs, seems to be 180° 

• e number of oscillations for this case, required 
,? r ™ e ' siufted “alysis is approximately twice that of 
S ° passa S e analysis, resulting in roughly equal CPU 


requirements. 

™^,°^ tbe above imputations were carried out on 
CKAY C-90 computer at the NAS facility of NASA. Fo 
the time-shifted analysis approximately 15 minutes wer 
required per oscillation cycle per passage. A total of 7, 
Mw of memory was required for in-core storage of bound 
ary condition data, whereas only 35 Mw were required i 
the ssds were used. 


CONCLUSIONS 

A time-shifted boundary condition has been success 
fully implemented into the TURBO- AE, an Euler/Navier 
Stokes based aeroelastic solver. Although, the direct ston 
method has been implemented, the large storage require- 
ments are not needed because the data is written and 
read from the solid-state input-output devices on Cray 
computer. The solver has been verified by applying it 
to a flat plate helical fan geometry. The results indicate 
that the time-shifted boundary conditions have been im- 
plemented satisfactorily and provide solutions that are in 
good agreement with other published results and the mul- 
tiple passage analysis. 

The analysis showed that significant computational 
tune can be saved using this procedure over the method 
crfstaclring the blade passages for analyzing the non-zero 
mPA - Even ttl0 ugfl the number of oscillation cycles re- 
quired for convergence are higher for the time-shifted 
boundapr conditions, the overall computational cost is re- 
duced significantly, since only one passage is used in the 
analysis. It was also found that for the present geome- 
try, the number of oscillations required for convergence 
were not strongly coupled to IBPA of analysis and for the 
cases analyzed, the number of oscillations required were 
of the order of seven to ten. However, one must be cau- 

tioned that this conclusion may not hold strictly for other 
geometries. 
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Figure 2: Unsteady pressure difference variation with 
chord at mid-span for 0° IBPA pitching oscillations 



Figure 1: Effect of time steps on convergence of work-per- 
cycle 


Figure 3: Unsteady pressure difference variation with 
chord at mid-span for 180° IBPA pitching oscillations 
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figure 4: Unsteady pressure difference variation with 
h.ord at mid-span for 0° IBPA plunging oscillations 



Figure 5: Unsteady pressure difference variation with 
chord at mid-span for 180° IBPA plunging oscillations 
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Figure 6: Comparison of work-per-cycle convergence for 
time-shifted analysis with multiple passage analysis 



Figure 7: Comparison of unsteady pressure diff erence 
variation at mid-span for 180° IBPA pit chin g os cill ations 






Figure 8: Comparison of unsteady pressure difference 
variation at mid-span for 180° IB PA plunging oscillations 


Figure 10: Comparison of unsteady pressure difference 
variation at mid-span for 90° IBPA pitching oscillations 




Figure 9: Comparison of unsteady pressure difference 
variation at mid-span for -90° IBPA pitching oscillations 


Figure 11: Comparison of work-per-cycle convergence his- 
tory for various IBPA pitching oscillations 
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Abstract 

This paper presents representative results from an 
aeroelastic code (TURBO-AE) based on an 
Euler / Navier-Stokes unsteady aerodynamic code 
(TURBO). Unsteady pressure, lift, and moment 
distributions are presented for a helical fan test 
configuration which is used to verify the code by 
comparison to two-dimensional linear potential (flat 
plate) theory. The results are for pitching and 
plunging motions over a range of phase angles. Good 
agreement with linear theory is seen for all phase 
angles except those near acoustic resonances. The 
agreement is better for pitching motions than for 
plunging motions. The reason for this difference is 
not understood at present. Numerical checks have 
been performed to ensure that solutions are 
independent of time step, converged to periodicity, 
and linearly dependent on amplitude of blade motion. 
The paper concludes with an evaluation of the 
current state of development of the TURBO-AE code 
and presents some plans for further development and 
validation of the TURBO-AE code. 
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Introduction 

There is an ongoing effort to develop technologies to 
increase the fuel efficiency of commercial aircraft 
engines, improve the safety of engine operation, 
reduce the emissions, and reduce engine noise. With 
the development of new designs of ducted fans, 
compressors, and turbines to achieve these goals, a 
basic aeroelastic requirement is that there should be 
no flutter or high resonant blade stresses in the 
operating regime. In order to verify the aeroelastic 
soundness of the design, an accurate prediction of the 
unsteady aerodynamics and structural dynamics of 
the propulsion component is required. The complex 
geometry, the presence of shock waves and flow 
separation makes the modeling of the unsteady 
aerodynamics a difficult task. The advanced blade 
geometry, new blade materials and new blade 
attachment concepts make the modeling of the 
structural dynamics a difficult problem. 

Computational aeroelastic modeling of fans, 
compressors, and turbines requires many simplifying 
assumptions. For instance, flutter calculations are 
typically carried out assuming that the blade row is 
isolated. This simplifies the structural dynamics 
formulation and the unsteady aerodynamic 
calculations considerably. 

For an isolated blade row flutter calculation, the 
modeling of the unsteady aerodynamics is the biggest 
challenge. Many simplifying assumptions are made 
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in the modeling of the unsteady aerodynamics. In the 
past, panel methods based on linear compressible 
small-disturbance potential theory have been used to 
model the unsteady aerodynamics and aeroelasticity 
of fans in subsonic flow; see for example [1,2]. The 
major limitations of this type of analysis are the 
neglect of transonic, vortical, and viscous flow effects 
in the model. These inherent limitations in the model 
preclude its use in a majority of practical 
applications. A full potential unsteady aerodynamic 
analysis has been used with a modal structural 
dynamics method to model the aeroelastic behavior of 
fan blades [3,4]. Although the full potential 
aerodynamic formulation is able to model transonic 
effects (limited to weak shocks), the vortical and 
viscous effects are still neglected. For example, the 
blade tip vortex, or a leading-edge vortex is not 
modeled. Recently, researchers [5-10] have also 
developed inviscid and viscous unsteady aerodynamic 
analyses for vibrating blades. 

For aeroelastic problems in which, viscous effects play 
an important role (such as flutter with flow 
separation, or stall flutter, and flutter in the presence 
of shock and boundary-layer interaction), a more 
advanced aeroelastic computational capability is 
required. The authors of this paper have earlier 
presented [11] some results from the TURBO-AE 
aeroelastic code. Initial calculations were restricted 
to in-phase (zero phase angle) blade motions and 
inviscid flow. In a later paper [12], results were 
presented for zero and non-zero phase angle motions 
and viscous flow. In these calculations, multiple 
blade passages were modeled for non-zero phase 
angle motions. Most recently [13], results have been 
presented using a single blade passage with phase- 
lag periodic boundary conditions to model arbitrary 
phase angle motions. 

This paper presents unsteady pressure, lift, and 
moment distributions due to blade vibration over a 
range of phase angles for verification of the TURBO- 
AE aeroelastic code. For non-zero phase angle 
motions, phase-lag periodic boundary conditions are 
used. The configuration selected is a helical fan. The 
geometry and flow conditions are chosen to minimize 
non-linear and three-dimensional effects since the 
intent is to verify the code by comparison with two- 
dimensional linear potential (flat plate) theory. 


Aeroelastic Code - TURBO-AE 

This section briefly describes the aeroelastic code 
(TURBO-AE); previous publications [11-13] provide 


additional details. The TURBO-AE code is based on 
an unsteady aerodynamic Euler / Navier-Stokes code 
(TURBO), developed separately [14,15]. The TURBO 
code provides all the unsteady aerodynamics to the 
TURBO-AE code. 

The TURBO code was originally developed [14] as an 
inviscid flow solver for modeling the flow through 
turbomachinery blade rows. Additional developments 
were made [15] to incorporate viscous effects into the 
model. This Reynolds-averaged Navier Stokes 
unsteady aerodynamic code is based on a finite 
volume scheme. Flux vector splitting is used to 
evaluate the flux Jacobians on the left band side of 
the governing equations [14] and Roe s flux difference 
splitting is used to form a higher-order TVD (Total 
Variation Diminishing) scheme to evaluate the fluxes 
on the right hand side. Newton sub -iterations are 
used at each time step to maintain higher accuracy. 
Symmetric Gauss-Seidel iterations are applied to the 
discretized equations. A Baldwin-Lomax algebraic 
turbulence model is used in the code. 

The TURBO-AE code assumes a normal mode 
representation of the structural dynamics of the 
blade. A work-per-cycle method is used to determine 
aeroelastic stability (flutter). Using this method, the 
motion of the blade is prescribed to be a harmonic 
vibration in a specified in-vacuum normal mode with 
a specified frequency (typically the natural 
frequency). The work done on the vibrating blade by 
aerodynamic forces during a cycle of vibration is 
calculated. If work is being done on the blade by the 
aerodynamic forces at the end of a vibration cycle, the 
blade is dynamically unstable, since it will result in 
extraction of energy from the flow, leading to an 
increase in amplitude of oscillation of the blade. 

The inlet/exit boundary conditions used in this code 
are described in [16-18], For cases in which the blade 
motions are not in-phase, phase-lag periodic 
boundary conditions based on the direct store method 
are used. 


Results 

In this section, results are presented which serve to 
verify the TURBO-AE code. The test configuration 
selected is a helical fan [16]. This configuration 
consists of a rotor with twisted flat plate blades 
enclosed in a cylindrical duct with no tip gap. This 
configuration was developed by researchers [16] to 
provide a relatively simple test case for comparison 
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with : T »vo-dimensionai analyses. The geometry' is such 
that three-dimensionality of the flow is minimized. 

The parameters of this three-dimensional 
configuration are such that the mid-span location 
corresponds to a flat plate cascade with a stagger 
angle of 45 deg. and unit gap-to-chord ratio operating 
in a uniform mean flow at a Mach number of 0.7 
parallel to the blades. The rotor has 24 blades with a 
hub/tip ratio of 0.8. The inlet flow (axial) Mach 
number used in this calculation is 0.495, which 
results in a relative Mach number of approximately 
0.7 at the mid-span section. The results presented 
are for inviscid runs of the TURBO-AE code. 

The grid used for the calculations is 141x11x41 in 
one blade passage. On each blade surface, 81 points 
are located in the chordwise direction and 11 points 
in the spanwise direction. The inlet and exit 
boundaries are located at an axial distance of 
approximately 0.7 chord lengths from the blade 
leading and trailing edges. To begin, a steady 
solution is obtained for this configuration. The steady 
flowfield consists of uniform flow at each radial 
location. 

Aeroelastic calculations are performed starting from 
the steady solution. Calculations have been 
performed for harmonic blade vibration in plunging 
and pitching modes, separately. The pitching is about 
the mid-chord. The prescribed mode shapes are such 
that the amplitude of vibration does not vary along 
the span. This choice of mode shapes is meant to 
reduces the three-dimensionality of the unsteady 
flowfield for ease of comparison with two-dimensional 
analyses. 

The vibration frequency is selected so that the non- 
dimensional reduced frequency based on blade chord 
is 1.0 at the mid-span. A study was performed to 
determine the sensitivity of numerical results to the 
number of time steps used in each cycle of blade 
vibration. Calculations were done with 100, 200, and 
300 time steps per cycle of vibration for 0 deg. phase 
angle plunging motion. The time step was varied so 
as to keep the vibration time period (or frequency) 
fixed. Figure 1 shows the work-per-cycle from this 
study. As the flowfield reaches periodicity, it can be 
seen that the results are nearly identical for 200 and 
300 time steps per cycle. These results differ slightly 
from the results for 100 time steps per cycle. Figure 
2 shows the unsteady pressure difference for the 
same three numbers of time steps per cycle. The 
results for 200 and 300 time steps per cycle are 


indistinguishable. Based on such calculations, it was 
determined that 200 time steps per cycle provided 
adequate temporal resolution for the selected 
vibration frequency. All results presented here have 
been obtained using 200 time steps per cycle . 

The non-dimensional time step used in the 
calculations (with 200 time steps per cycle) is 0.045, 
which results in a maximum CFL number of 60.5. 
The amplitude of blade vibrations in the calculation 
is a pitching amplitude of 0.2 deg. or a plunging 
amplitude of 0.1% chord. In all cases, calculations 
were continued for a number of cycles of blade 
vibration to allow the flowfield to become periodic. 
Initial calculations with phase angles of 0, 45, 90, 
135, 180, 225, 270, and 315 deg. were continued for 
15 cycles of blade vibration to ensure periodicity. 
Later calculations with intermediate phase angles 
(22.5, 67.5, ... , and 337.5 deg.) were continued only 
for 10 cycles of blade vibration due to insufficient 
computational resources. In an earlier study [13], it 
was shown that, for the various phase angles studied, 
the flowfield became periodic after about 7-10 cycles 
of blade vibration. Hence, the 10 or 15 cycles used in 
the present work were considered adequate to reach 
periodicity. 

Figure 3 shows the unsteady moment about mid- 
chord (in complex form) for pitching blade motion 
about the mid-chord. These results are from the mid- 
span location and were calculated using the first 
harmonic of the unsteady blade surface pressure 
difference. Semi-analytical results from two- 
dimensional linear potential (flat plate) theory [19] 
are included for comparison. 

The overall level of agreement between TURBO-AE 
results and linear theory is very good, with 
exceptions to be discussed in the following paragraph. 
For subsonic flows and small amplitude of blade 
motions, it is expected that there will be no 
significant difference between the Euler and linear 
potential results. Hence, the observed agreement is 
not surprising and provides a basic verification of the 
TURBO-AE code. It may be noted that the 
parameters of the present configuration were 
selected [16] to allow exactly this type of a 
verification by comparison to two-dimensional 
analyses. 

In Figure 3, some deviation from linear theory is seen 
in the results for phase angles of 112.5 and 135 deg., 
and to a lesser extent for phase angles of 157.5 and 
315 deg. All these phase angles fall near conditions of 
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acoustic resonance (or cut-off conditions) in the 
corresponding two-dimensional flat plate cascade. 
The acoustic resonances occur at phase angles of 
107.3 and 330.6 deg.; these values are marked on the 
phase angle axis of Figure 3 for reference. The phase 
angles between these resonances are associated with 
sub-resonant [20] (cut-off) conditions in which all 
disturbances attenuate away from the cascade. No 
disturbances propagate in the upstream or 
downstream directions under sub-resonant 
conditions. The phase angles between 0 and 107.3 
deg. and between 330.6 and 360 deg. are associated 
with super-resonant (cut-on) conditions in which at 
least one disturbance propagates in either the far 
upstream or downstream direction. 

The significance of the sub-resonant and super- 
resonant conditions to computational aeroelasticity 
can be explained as follows. Since the typical 
computational domain does not extend very far from 
the blade row or cascade, the inlet/exit boundary 
conditions must minimize (or eliminate) the 
reflection of disturbances generated by the vibration 
of the blades. For sub-resonant conditions, it may be 
possible to reduce the reflected disturbances by 
moving the boundary farther away from the blade 
row. This is not possible for super-resonant 
conditions. From Figure 3, it can be seen that the 
results from TURBO- AE agree well with linear 
theory for both sub-resonant and super-resonant 
conditions. It may be also recalled that the 
computational inlet/exit boundaries are located quite 
near (0.7 axial chord lengths from leading/trailing 
edges) the blade row in the present calculations. 

Figure 4 shows the unsteady lift (in complex form) for 
plunging blade motion. As noted for the pitching 
results, these results are also from the mid-span 
location and were also calculated using the first 
harmonic of the unsteady blade surface pressure 
difference. Results from linear potential theory are 
included in Figure 4 for comparison. The overall level 
of agreement with linear theory is good, but not as 
good as that for pitching motion (Figure 3). The 
source of such a difference between the plunging and 
pitching results is not understood. However, such 
differences in agreement have been noted by other 
researchers [16,17] for a different configuration. In 
addition, deviations are observed close to the acoustic 
resonances, as for pitching. 

Figure 5 shows the unsteady blade surface pressure 
difference (first harmonic) at the mid-span location 
for pitching blade motion about the mid-chord. 


Results are presented for phase angles values 
between 0 and 360 deg. in steps of 22.5 deg. In each 
case, the linear theory results are included for 
comparison. In most cases, the agreement with linear 
theory is very good. The exceptions occur at phase 
angles near acoustic resonance conditions, as noted 
earlier in the description of the unsteady moment 
(Figure 3). It is worth noting that, in this case, the 
integrated results in Figure 3 accurately represent 
the level of agreement with linear theory, without 
obscuring any differences in the details of the 
pressure distributions. 

Figure 6 shows the unsteady blade surface pressure 
difference (in complex form) for plunging blade 
motion. The level of agreement with linear theory is 
not as good as for pitching, as reflected in the 
unsteady lift (Figure 4). The most serious deviations 
from linear theory are restricted to the phase angles 
near conditions of acoustic resonance. 

Some of the results for plunging motion (Figure 6) 
show an irregular (unsmooth) variation in the 
unsteady pressure distribution which is not seen in 
any of the results for pitching motion (Figure 5). This 
uneven variation can be seen in the plunging results 
in Figures 6b, 6d, 6f, 6h, 6j, 61, 6n, and 6p for phase 
angles of 22.5, 67.5, 112.5, ... , and 337.5 deg. One 
common characteristic of these results is that these 
were all generated on a workstation and may 
therefore suffer from some precision-related 
numerical problem. However, it is surprising to note 
that the corresponding results for pitching motion 
(also computed on a workstation) are quite smooth 
and do not show such unevenness. A re-calculation of 
selected plunging results on a super-computer does 
indeed eliminate the unevenness in pressure 
variation, but the pressure distributions remain 
substantially unchanged from those presented in 
Figure 6. 

Note that all the TURBO-AE results presented are 
the first harmonic components of the unsteady 
variations. The higher harmonics are extremely 
small for these calculations, indicating the linearity 
of the unsteady flow. Previous results [12] had shown 
a nonlinear dependence on amplitude for certain 
cases for pitching amplitudes of blade vibration of 2 
deg., but not at the 0.2 deg. amplitude used in the 
present calculations. 

To investigate the effect of some numerical 
parameters on the results for phase angle of 112.5 
deg. (where the maximum deviation from linear 
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theory is observed], the rollowing calculations were 
done. The number of time steps per cycle was doubled 
from 200 to 400, with a corresponding halving of the 
time step. The unsteady pressure results showed no 
changes within plotting accuracy, indicating 
adequate temporal resolution. Similarly, the number 
of cycles of oscillation was doubled from 10 to 20 to 
examine possible lack of periodicity. No change in the 
unsteady pressure results was observed within 
plotting accuracy. The deviations in the regions of 
acoustic resonances may possibly be reduced by the 
use of finer grids. But, such a grid refinement study 
has not yet been performed. 


Concluding Remarks 

An aeroelastic analysis code named TURBO-AE has 
been developed and is being verified and validated. 
The starting point for the development was an 
Euler / Navier-Stokes unsteady aerodynamic code 
named TURBO. Some verification has been done by 
running the code for a helical fan test configuration. 
Results have been presented for pitching and 
plunging blade motions over a range of phase angles. 
The results compare well with results from a linear 
potential analysis. This agreement is expected for 
subsonic flows for which the calculations were made 
and for the relatively small amplitudes of blade 
motion. 

The agreement is not as good for plunging motion as 
for pitching motion. The reason for this difference is 
not understood at present. Also, deviations are 
observed for values of phase angles near acoustic 
resonance conditions. The solutions are shown to be 
independent of the time step, converged to 
periodicity, and linearly dependent on amplitude of 
blade motion. This test case provides a basic 
verification of the TURBO-AE code. It also shows the 
need to perform a grid refinement study as a possible 
way to resolve the deviations from linear theory near 
acoustic resonance conditions and for plunging 
motion. For plunging motion, some results are 
affected by precision-related numerical problems, as 
seen from uneven pressure distributions. But, the 
elimination of these precision problems does not 
change the pressure distributions substantially, 
apart from making the variations smooth. 

It is necessary to further verify the TURBO-AE using 
different standard test configurations to compare 
with experimental data and other code predictions. 


This is being done in collaboration with other 
researchers. Also, it is necessary that the TURBO-AE 
code be exercised to evaluate its ability to analyze 
and predict flutter for conditions in which viscous 
effects are significant. This work is also currently in 
progress. 
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Figure 3: Unsteady moment for pitching motion. 
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Figure 4: Unsteady lift for plunging motion. 
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(a) 0 deg. pitching 


(e) 90 deg. pitching 




(f) 112.5 deg. pitching 



(c) 45 deg. pitching 



(g) 135 deg. pitching 




(h) 157.5 deg. pitching 


Figure 5. Unsteady pressure difference (first harmonic) for pitching motion. 
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Unsteady pressure difference Unsteady pressure difference Unsteady pressure difference Unsteady pressure difference 




(i) 180 deg. pitching 


(m) 270 deg. pitching 




(j) 202.5 deg. pitching (n) 292.5 deg. pitching 




(k) 225 deg. pitching (o) 315 deg. pitching 




(1) 247.5 deg. pitching (p) 337.5 deg. pitching 


Figure 5 (continued): Unsteady pressure difference (first harmonic) for pitching motion. 
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(b) 22.5 deg. plunging 


(f) 112.5 deg. plunging 




(c) 45 deg. plunging 


(g) 135 deg. plunging 




(d) 67.5 deg. plunging (h) 157.5 deg. plunging 


Figure 6: Unsteady pressure difference (first harmonic) for plunging motion. 
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Unsteady pressure difference Unsteady pressure difference Unsteady pressure difference Unsteady pressure diffe 



(i) 180 deg. plunging 


(m) 270 deg. plunging 



0‘) 202.5 deg. plunging (n) 292.5 deg. plunging 



(k) 225 deg. plunging (o) 315 deg. plunging 



(1) 247.5 deg. plunging (p) 337.5 deg. plunging 

Figure 6 (continued): Unsteady pressure difference (first harmonic) for plunging motion. 
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